{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Protein-ligand Docking tutorial using BioExcel Building Blocks (biobb)\n",
    "### -- *Fpocket Version* --\n",
    "\n",
    "***\n",
    "This tutorial aims to illustrate the process of **protein-ligand docking**, step by step, using the **BioExcel Building Blocks library (biobb)**. The particular example used is the **Mitogen-activated protein kinase 14** (p38-α) protein (PDB code [3HEC](https://www.rcsb.org/structure/3HEC)), a well-known **Protein Kinase enzyme**, \n",
    " in complex with the FDA-approved **Imatinib**, (PDB Ligand code [STI](https://www.rcsb.org/ligand/STI), DrugBank Ligand Code [DB00619](https://go.drugbank.com/drugs/DB00619)), a small molecule **kinase inhibitor** used to treat certain types of **cancer**. \n",
    " \n",
    "The tutorial will guide you through the process of identifying the **active site cavity** (pocket) without previous knowledge, and the final prediction of the **protein-ligand complex**. \n",
    "\n",
    "Please note that **docking algorithms**, and in particular, **AutoDock Vina** program used in this tutorial, are **non-deterministic**. That means that results obtained when running the workflow **could be diferent** from the ones we obtained during the writing of this tutorial (see [AutoDock Vina manual](http://vina.scripps.edu/manual.html)). We invite you to try the docking process several times to verify this behaviour. \n",
    "***\n",
    "\n",
    "\n",
    "<div style=\"background:#b5e0dd; padding: 15px;\"><strong>Important:</strong> it is recommended to execute this tutorial step by step (not as a single workflow execution, <strong><em>Run All</em></strong> mode), as it has interactive selections.</div>\n",
    "\n",
    "## Settings\n",
    "\n",
    "### Biobb modules used\n",
    "\n",
    " - [biobb_io](https://github.com/bioexcel/biobb_io): Tools to fetch biomolecular data from public databases.\n",
    " - [biobb_structure_utils](https://github.com/bioexcel/biobb_structure_utils): Tools to modify or extract information from a PDB structure file.\n",
    " - [biobb_chemistry](https://github.com/bioexcel/biobb_chemistry): Tools to perform chemoinformatics processes.\n",
    " - [biobb_vs](https://github.com/bioexcel/biobb_vs): Tools to perform virtual screening studies.\n",
    " \n",
    "### Auxiliar libraries used\n",
    "\n",
    " - [nb_conda_kernels](https://github.com/Anaconda-Platform/nb_conda_kernels): Enables a Jupyter Notebook or JupyterLab application in one conda environment to access kernels for Python, R, and other languages found in other environments.\n",
    " - [nglview](http://nglviewer.org/#nglview): Jupyter/IPython widget to interactively view molecular structures and trajectories in notebooks.\n",
    " - [ipywidgets](https://github.com/jupyter-widgets/ipywidgets): Interactive HTML widgets for Jupyter notebooks and the IPython kernel.\n",
    "\n",
    "### Conda Installation and Launch\n",
    "\n",
    "```console\n",
    "git clone https://github.com/bioexcel/biobb_wf_virtual-screening.git\n",
    "cd biobb_wf_virtual-screening\n",
    "conda env create -f conda_env/environment.yml\n",
    "conda activate biobb_VS_tutorial\n",
    "jupyter-nbextension enable --py --user widgetsnbextension\n",
    "jupyter-nbextension enable --py --user nglview\n",
    "jupyter-notebook biobb_wf_virtual-screening/notebooks/fpocket/wf_vs_fpocket.ipynb\n",
    "  ``` \n",
    "\n",
    "***\n",
    "## Pipeline steps\n",
    " 1. [Input Parameters](#input)\n",
    " 2. [Fetching PDB Structure](#fetch)\n",
    " 3. [Extract Protein Structure](#extractProtein)\n",
    " 4. [Computing Protein Cavities (fpocket)](#fpocket)\n",
    " 5. [Filtering Protein Cavities (fpocket output)](#fpocketFilter)\n",
    " 6. [Extract Pocket Cavity ](#fpocketSelect)\n",
    " 7. [Generating Cavity Box ](#cavityBox)\n",
    " 8. [Downloading Small Molecule](#downloadSmallMolecule)\n",
    " 9. [Converting Small Molecule](#sdf2pdb)\n",
    " 10. [Preparing Small Molecule (ligand) for Docking](#ligand_pdb2pdbqt)\n",
    " 11. [Preparing Target Protein for Docking](#protein_pdb2pdbqt)\n",
    " 12. [Running the Docking](#docking)\n",
    " 13. [Extract a Docking Pose](#extractPose)\n",
    " 14. [Converting Ligand Pose to PDB format](#pdbqt2pdb)\n",
    " 15. [Superposing Ligand Pose to the Target Protein Structure](#catPdb)\n",
    " 16. [Comparing final result with experimental structure](#viewFinal)\n",
    " 17. [Questions & Comments](#questions)\n",
    " \n",
    "***\n",
    "<img src=\"https://bioexcel.eu/wp-content/uploads/2019/04/Bioexcell_logo_1080px_transp.png\" alt=\"Bioexcel2 logo\"\n",
    "\ttitle=\"Bioexcel2 logo\" width=\"400\" />\n",
    "***\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"input\"></a>\n",
    "## Input parameters\n",
    "**Input parameters** needed:\n",
    "\n",
    " - **pdb_code**: PDB code of the experimental complex structure (if exists).<br>\n",
    "In this particular example, the **p38α** structure in complex with the **Imatinib drug** was experimentally solved and deposited in the **PDB database** under the **3HEC** PDB code. The protein structure from this PDB file will be used as a **target protein** for the **docking process**, after stripping the **small molecule**. An **APO structure**, or any other structure from the **p38α** [cluster 100](https://www.rcsb.org/search?request=%7B%22query%22%3A%7B%22type%22%3A%22terminal%22%2C%22service%22%3A%22sequence%22%2C%22parameters%22%3A%7B%22target%22%3A%22pdb_protein_sequence%22%2C%22value%22%3A%22RPTFYRQELNKTIWEVPERYQNLSPVGSGAYGSVCAAFDTKTGLRVAVKKLSRPFQSIIHAKRTYRELRLLKHMKHENVIGLLDVFTPARSLEEFNDVYLVTHLMGADLNNIVKCQKLTDDHVQFLIYQILRGLKYIHSADIIHRDLKPSNLAVNEDCELKILDFGLARHTDDEMTGYVATRWYRAPEIMLNWMHYNQTVDIWSVGCIMAELLTGRTLFPGTDHIDQLKLILRLVGTPGAELLKKISSESARNYIQSLTQMPKMNFANVFIGANPLAVDLLEKMLVLDSDKRITAAQALAHAYFAQYHDPDDEPVADPYDQSFESRDLLIDEWKSLTYDEVISFVPPP%22%2C%22identity_cutoff%22%3A1%2C%22evalue_cutoff%22%3A0.1%7D%2C%22node_id%22%3A0%7D%2C%22return_type%22%3A%22polymer_entity%22%2C%22request_options%22%3A%7B%22pager%22%3A%7B%22start%22%3A0%2C%22rows%22%3A25%7D%2C%22scoring_strategy%22%3A%22combined%22%2C%22sort%22%3A%5B%7B%22sort_by%22%3A%22score%22%2C%22direction%22%3A%22desc%22%7D%5D%7D%2C%22request_info%22%3A%7B%22src%22%3A%22ui%22%2C%22query_id%22%3A%22bea5861f8b38a9e25a3e626b39d6bcbf%22%7D%7D) (sharing a 100% of sequence similarity with the **p38α** structure) could also be used as a **target protein**. This structure of the **protein-ligand complex** will be also used in the last step of the tutorial to check **how close** the resulting **docking pose** is from the known **experimental structure**. \n",
    " -----\n",
    " - **ligandCode**: Ligand PDB code (3-letter code) for the small molecule (e.g. STI).<br>\n",
    "In this particular example, the small molecule chosen for the tutorial is the FDA-approved drug **Imatinib** (PDB Code STI), a type of cancer growth blocker, used in [diferent types of leukemia](https://go.drugbank.com/drugs/DB00619).\n",
    " -----\n",
    " - **pockets_dir**: Name of a folder to write temporary files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nglview\n",
    "import ipywidgets\n",
    "\n",
    "pdb_code = \"3HEC\"         # P38 + Imatinib\n",
    "\n",
    "ligand_code = \"STI\"       # Imatinib\n",
    "\n",
    "pockets_dir = \"pockets\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"fetch\"></a>\n",
    "***\n",
    "## Fetching PDB structure\n",
    "Downloading **PDB structure** with the **protein molecule** from the PDBe database.<br>\n",
    "Alternatively, a **PDB file** can be used as starting structure. <br>\n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [Pdb](https://biobb-io.readthedocs.io/en/latest/api.html#module-api.pdb) from **biobb_io.api.pdb**\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from biobb_io.api.pdb import pdb\n",
    "\n",
    "download_pdb = \"download.pdb\"\n",
    "prop = {\n",
    "  \"pdb_code\": pdb_code,\n",
    "  \"filter\": [\"ATOM\", \"HETATM\"]\n",
    "}\n",
    "\n",
    "pdb(output_pdb_path=download_pdb,\n",
    "    properties=prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"vis3D\"></a>\n",
    "### Visualizing 3D structure\n",
    "Visualizing the downloaded/given **PDB structure** using **NGL**.<br><br>\n",
    "Note (and try to identify) the **Imatinib small molecule (STI)** and the **detergent (β-octyl glucoside) (BOG)** used in the experimental reservoir solution to obtain the crystal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "view = nglview.show_structure_file(download_pdb, default=True)\n",
    "view.center()\n",
    "view._remote_call('setSize', target='Widget', args=['','600px'])\n",
    "\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='ngl1.png'></img>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"extractProtein\"></a>\n",
    "***\n",
    "## Extract Protein Structure\n",
    "Extract **protein structure** from the **downloaded PDB file**. Removing **any extra molecule** (ligands, ions, water molecules). <br><br>\n",
    "The **protein structure** will be used as a **target** in the **protein-ligand docking process**. \n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [extract_molecule](https://biobb-structure-utils.readthedocs.io/en/latest/utils.html#module-utils.extract_molecule) from **biobb_structure_utils.utils.extract_molecule**\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_structure_utils.utils.extract_molecule import extract_molecule\n",
    "\n",
    "pdb_protein = \"pdb_protein.pdb\"\n",
    "\n",
    "extract_molecule(input_structure_path=download_pdb,\n",
    "             output_molecule_path = pdb_protein)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"vis3D\"></a>\n",
    "### Visualizing 3D structure\n",
    "Visualizing the downloaded/given **PDB structure** using **NGL**.<br><br>\n",
    "Note that the **small molecules** included in the original structure are now gone. The new structure only contains the **protein molecule**, which will be used as a **target** for the **protein-ligand docking**. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "view = nglview.show_structure_file(pdb_protein, default=False)\n",
    "view.add_representation(repr_type='cartoon', \n",
    "                        selection='not het',\n",
    "                       colorScheme = 'atomindex')\n",
    "view.center()\n",
    "view._remote_call('setSize', target='Widget', args=['','600px'])\n",
    "\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='ngl2.png'></img>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"fpocket\"></a>\n",
    "***\n",
    "## Computing Protein Cavities (fpocket)\n",
    "Computing the **protein cavities** (pockets) using the well-known [**fpocket**](https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-10-168) tool.<br>\n",
    "These **cavities** will be then used in the **docking procedure** to try to find the **best region of the protein surface** where the small molecule can **bind**. <br><br>\n",
    "Although in this particular example we already know the **binding site** region, as we started from a **protein-ligand complex** structure where the ligand was located in the same **binding site** as **Imatinib** is binding, this is not always the case. In the cases where we do not know these regions, **fpocket** will help us identifying the possible **binding sites** of our **target protein**.<br>\n",
    "\n",
    "**fpocket** input parameters, such as **minimum** and **maximum radius** (in Angstroms) the alpha spheres might have in a **binding pocket** can be adjusted (min_radius, max_radius) . Parameters used in this particular example are 3Å for the **minimum radius** and 6Å for the **maximum radius**. The **minimum number of alpha spheres** a pocket must contain in order to figure in the results is also adjusted to 35. See the [fpocket manual](http://fpocket.sourceforge.net/manual_fpocket2.pdf) for more information.<br> \n",
    "<br>\n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [fpocket_run](https://biobb-vs.readthedocs.io/en/latest/fpocket.html#module-fpocket.fpocket_run) from **biobb_vs.fpocket.fpocket_run**\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from biobb_vs.fpocket.fpocket_run import fpocket_run\n",
    "\n",
    "fpocket_all_pockets = \"fpocket_all_pockets.zip\"\n",
    "fpocket_summary = \"fpocket_summary.json\"\n",
    "prop = {\n",
    "    \"min_radius\": 3,\n",
    "    \"max_radius\": 6,\n",
    "    \"num_spheres\": 35\n",
    "}\n",
    "\n",
    "fpocket_run(input_pdb_path=pdb_protein,\n",
    "        output_pockets_zip = fpocket_all_pockets,\n",
    "        output_summary=fpocket_summary,\n",
    "        properties=prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"checkJson\"></a>\n",
    "### Checking fpocket output (json)\n",
    "Checking the **fpocket** output from the **json file**. Every **pocket** has a separated entry in the json output, with information such as: **score, druggability score, volume, hydrophobicity, polarity or flexibility**. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import json\n",
    "\n",
    "with open(fpocket_summary, 'r') as json_file:\n",
    "    data = json.load(json_file)\n",
    "    print(json.dumps(data, indent=4))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"fpocketFilter\"></a>\n",
    "***\n",
    "## Filtering Protein Cavities (fpocket output)\n",
    "Filtering the **protein cavities** (pockets) identified by **fpocket**.<br>\n",
    "In this particular example, the biggest **cavities**, with a **volume** between 800 and 2000 ($Å^{3}$), big enough volume to fit the input small molecule, are selected. <br>\n",
    "\n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [fpocket_filter](https://biobb-vs.readthedocs.io/en/latest/fpocket.html#module-fpocket.fpocket_filter) from **biobb_vs.fpocket.fpocket_filter**\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from biobb_vs.fpocket.fpocket_filter import fpocket_filter\n",
    "\n",
    "fpocket_filter_pockets = \"fpocket_filter_pockets.zip\"\n",
    "prop = {\n",
    "    \"volume\": [800, 2000]\n",
    "}\n",
    "\n",
    "fpocket_filter(input_pockets_zip=fpocket_all_pockets,\n",
    "                input_summary = fpocket_summary,\n",
    "                output_filter_pockets_zip=fpocket_filter_pockets,\n",
    "                properties=prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"extractPockets\"></a>\n",
    "### Extract selected pockets (cavities)\n",
    "Extract the selected **pockets** (cavities) from the filtered list (zip file, fpocket_filter_pockets).<br>\n",
    "Writing the information in the ***pockets_dir*** folder.<br>\n",
    "Also saving the list of **PDB files** (protein residues forming the pocket) and **PQR files** (cavity, pocket), to be used in following **visualization step**. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import shutil\n",
    "\n",
    "from pathlib import Path, PurePath\n",
    "import zipfile\n",
    "\n",
    "if Path(pockets_dir).exists(): shutil.rmtree(pockets_dir) \n",
    "os.mkdir(pockets_dir)\n",
    "\n",
    "with zipfile.ZipFile(fpocket_filter_pockets, 'r') as zip_ref:\n",
    "    zip_ref.extractall(pockets_dir)\n",
    "\n",
    "path_pockets = [str(i) for i in Path(pockets_dir).iterdir()]\n",
    "path_pockets_pdb = [str(i) for i in Path(pockets_dir).iterdir() if PurePath(i).suffix == '.pdb']\n",
    "path_pockets_pqr = [str(i) for i in Path(pockets_dir).iterdir() if PurePath(i).suffix == '.pqr']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"viewPockets\"></a>\n",
    "### Visualizing selected pockets (cavities)\n",
    "Visualizing the selected **pockets** (cavities) from the filtered list using **NGL viewer**.<br>\n",
    "\n",
    "**Protein residues** forming the **cavity** are represented in **random-colored surfaces**. **Pockets** are represented in a **blue-colored mesh**. Different **pockets** are identified with a floating **label**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import re\n",
    "import random\n",
    "\n",
    "# random colors for cavities\n",
    "r = lambda: random.randint(0,255)\n",
    "\n",
    "# load structure\n",
    "view = nglview.NGLWidget()\n",
    "c = view.add_component(pdb_protein)\n",
    "\n",
    "# load cavities (d) and pockets (p) and create pocketNames list\n",
    "c = {}\n",
    "p = {}\n",
    "pocketNames = []\n",
    "for pock in path_pockets:\n",
    "    g = re.findall('(?:pocket)(\\d+)(?:_\\w+)\\.(\\w+)', pock)\n",
    "    i = g[0][0]\n",
    "    suff = g[0][1]\n",
    "    if not [item for item in pocketNames if ('pocket' + i) in item]: pocketNames.append(('pocket' + i, int(i)))\n",
    "\n",
    "    if suff == 'pdb':\n",
    "        c[i] = view.add_component(pock)\n",
    "        c[i].clear()\n",
    "    else:\n",
    "        p[i] = view.add_component(pock)\n",
    "        p[i].clear()\n",
    "\n",
    "# sort pocket names\n",
    "pocketNames.sort(key=lambda tup: tup[1])\n",
    "        \n",
    "# representation for cavities\n",
    "for pock in path_pockets_pdb:\n",
    "    g = re.findall('(?:pocket)(\\d+)(?:_\\w+)\\.(\\w+)', pock)\n",
    "    i = g[0][0]\n",
    "    c[i].add_surface(color='#%02X%02X%02X' % (r(),r(),r()),  \n",
    "                     radius='1.5',\n",
    "                     lowResolution= True,\n",
    "                     # 0: low resolution \n",
    "                     smooth=1,\n",
    "                     useWorker= True,\n",
    "                     wrap= True)\n",
    "    \n",
    "# representation for pockets\n",
    "for pock in path_pockets_pqr:\n",
    "    g = re.findall('(?:pocket)(\\d+)(?:_\\w+)\\.(\\w+)', pock)\n",
    "    i = g[0][0]\n",
    "    p[i].add_surface( component=i, color='skyblue', surfaceType= 'av', contour=True )\n",
    "\n",
    "view.center()\n",
    "view._remote_call('setSize', target='Widget', args=['','600px'])\n",
    "view\n",
    "\n",
    "# show pocket labels\n",
    "code = \"\"\"\n",
    "var stage = this.stage;\n",
    "var view = this.stage.viewer;\n",
    "var clist_len = stage.compList.length;\n",
    "var i = 0;\n",
    "for(i = 0; i <= clist_len; i++){\n",
    "    if(stage.compList[i] != undefined && stage.compList[i].structure != undefined && stage.compList[i].object.name.indexOf('pqr') != -1) {        \n",
    "\n",
    "        var elm = document.createElement(\"div\");\n",
    "        elm.innerText = 'pocket' + stage.compList[i].object.name.match(/\\d+/g)\n",
    "        elm.style.color = \"black\";\n",
    "        elm.style.background = \"rgba(201, 149, 6, .8)\";\n",
    "        elm.style.padding = \"8px\";\n",
    "        \n",
    "        stage.compList[i].addAnnotation(stage.compList[i].structure.center, elm)\n",
    "    }\n",
    "}\n",
    "\"\"\"\n",
    "\n",
    "view._execute_js_code(code)\n",
    "\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='ngl3.png'></img>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"selectPockets\"></a>\n",
    "### Select pocket (cavity)\n",
    "Select a specific **pocket** (cavity) from the filtered list to be used in the **docking procedure**. <br>\n",
    "\n",
    "If **fpocket** has been able to identify the correct **binding site**, which we know from the original **protein-ligand structure**, it just needs to be selected. In this particular example, the pocket we are interested in is the **pocket number 6**. <br>\n",
    "\n",
    "Choose a **pocket** from the **DropDown list**:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mdsel = ipywidgets.Dropdown(\n",
    "    options=pocketNames,\n",
    "    description='Sel. pocket:',\n",
    "    disabled=False,\n",
    ")\n",
    "display(mdsel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"fpocketSelect\"></a>\n",
    "***\n",
    "## Extract Pocket Cavity \n",
    "Extract the selected **protein cavity** (pocket) from the **fpocket** results.<br>\n",
    "\n",
    "It will be used to generate the **docking box** in the **following step**.\n",
    "\n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [fpocket_select](https://biobb-vs.readthedocs.io/en/latest/fpocket.html#module-fpocket.fpocket_select) from **biobb_vs.fpocket.fpocket_select**\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_vs.fpocket.fpocket_select import fpocket_select\n",
    "\n",
    "fpocket_cavity = \"fpocket_cavity.pdb\"\n",
    "fpocket_pocket = \"fpocket_pocket.pqr\"\n",
    "prop = {\n",
    "    \"pocket\": mdsel.value\n",
    "}\n",
    "\n",
    "fpocket_select(input_pockets_zip=fpocket_filter_pockets,\n",
    "                output_pocket_pdb = fpocket_cavity,\n",
    "                output_pocket_pqr=fpocket_pocket,\n",
    "                properties=prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"cavityBox\"></a>\n",
    "***\n",
    "## Generating Cavity Box \n",
    "Generating a **box** surrounding the selected **protein cavity** (pocket), to be used in the **docking procedure**. The **box** is defining the region on the **surface** of the **protein target** where the **docking program** should explore a possible **ligand dock**.<br>\n",
    "An offset of **12 Angstroms** is used to generate a **big enough box** to fit the **small molecule** and its possible rotations.<br>\n",
    "\n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [box](https://biobb-vs.readthedocs.io/en/latest/utils.html#module-utils.box) from **biobb_vs.utils.box**\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_vs.utils.box import box\n",
    "\n",
    "output_box = \"box.pdb\"\n",
    "prop = {\n",
    "    \"offset\": 12,\n",
    "    \"box_coordinates\": True\n",
    "}\n",
    "\n",
    "box(input_pdb_path = fpocket_pocket,\n",
    "            output_pdb_path = output_box,\n",
    "            properties=prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"vis3D\"></a>\n",
    "### Visualizing binding site box in 3D structure\n",
    "Visualizing the **protein structure**, the **selected cavity**, and the **generated box**, all together using **NGL** viewer. Using the **original structure** with the **small ligand** inside (Imatinib, [STI](https://www.rcsb.org/ligand/STI)), to check that the **selected cavity** is placed in the **same region** as the **original ligand**. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "view = nglview.NGLWidget()\n",
    "s = view.add_component(download_pdb)\n",
    "b = view.add_component(output_box)\n",
    "p = view.add_component(fpocket_pocket)\n",
    "p.clear()\n",
    "\n",
    "atomPair = [\n",
    "    [ \"9999:Z.ZN1\", \"9999:Z.ZN2\" ],\n",
    "    [ \"9999:Z.ZN2\", \"9999:Z.ZN4\" ],\n",
    "    [ \"9999:Z.ZN4\", \"9999:Z.ZN3\" ],\n",
    "    [ \"9999:Z.ZN3\", \"9999:Z.ZN1\" ],\n",
    "    \n",
    "    [ \"9999:Z.ZN5\", \"9999:Z.ZN6\" ],\n",
    "    [ \"9999:Z.ZN6\", \"9999:Z.ZN8\" ],\n",
    "    [ \"9999:Z.ZN8\", \"9999:Z.ZN7\" ],\n",
    "    [ \"9999:Z.ZN7\", \"9999:Z.ZN5\" ],\n",
    "    \n",
    "    [ \"9999:Z.ZN1\", \"9999:Z.ZN5\" ],\n",
    "    [ \"9999:Z.ZN2\", \"9999:Z.ZN6\" ],\n",
    "    [ \"9999:Z.ZN3\", \"9999:Z.ZN7\" ],\n",
    "    [ \"9999:Z.ZN4\", \"9999:Z.ZN8\" ]\n",
    "]\n",
    "\n",
    "# structure\n",
    "s.add_representation(repr_type='cartoon', \n",
    "                        selection='not het',\n",
    "                        color='#cccccc',\n",
    "                       opacity=.2)\n",
    "# ligands box\n",
    "b.add_representation(repr_type='ball+stick',\n",
    "                     selection='9999',\n",
    "                     color='pink', \n",
    "                     aspectRatio = 8)\n",
    "# lines box\n",
    "b.add_representation(repr_type='distance', \n",
    "                     atomPair= atomPair,\n",
    "                     labelVisible=False,\n",
    "                     color= 'black')\n",
    "\n",
    "# pocket\n",
    "p.add_surface(component=mdsel.value, \n",
    "              color='skyblue', \n",
    "              surfaceType= 'av', \n",
    "              lowResolution= True,\n",
    "              # 0: low resolution \n",
    "              smooth=1,\n",
    "              contour=True,\n",
    "              opacity=0.4,\n",
    "              useWorker= True,\n",
    "              wrap= True )\n",
    "\n",
    "\n",
    "view.center()\n",
    "view._remote_call('setSize', target='Widget', args=['','600px'])\n",
    "\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='ngl4.png'></img>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"downloadSmallMolecule\"></a>\n",
    "***\n",
    "## Downloading Small Molecule \n",
    "Downloading the desired **small molecule** to be used in the **docking procedure**. <br>\n",
    "In this particular example, the small molecule of interest is the FDA-approved drug **Imatinib**, with PDB code **STI**.<br>\n",
    "\n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [ideal_sdf](https://biobb-io.readthedocs.io/en/latest/api.html#module-api.ideal_sdf) from **biobb_io.api.ideal_sdf**\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_io.api.ideal_sdf import ideal_sdf\n",
    "\n",
    "sdf_ideal = \"ideal.sdf\"\n",
    "prop = {\n",
    "  \"ligand_code\": ligand_code\n",
    "}\n",
    "\n",
    "ideal_sdf(output_sdf_path=sdf_ideal,\n",
    "    properties=prop)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"sdf2pdb\"></a>\n",
    "***\n",
    "## Converting Small Molecule \n",
    "Converting the desired **small molecule** to be used in the **docking procedure**, from **SDF** format to **PDB** format using the **OpenBabel chemoinformatics** tool. <br>\n",
    "\n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [babel_convert](https://biobb-chemistry.readthedocs.io/en/latest/babelm.html#module-babelm.babel_convert) from **biobb_chemistry.babelm.babel_convert**\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_chemistry.babelm.babel_convert import babel_convert\n",
    "\n",
    "ligand = \"ligand.pdb\"\n",
    "prop = {\n",
    "    \"input_format\": \"sdf\",\n",
    "    \"output_format\": \"pdb\",\n",
    "    \"obabel_path\": \"obabel\"\n",
    "}\n",
    "\n",
    "babel_convert(input_path = sdf_ideal,\n",
    "            output_path = ligand,\n",
    "            properties=prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"ligand_pdb2pdbqt\"></a>\n",
    "***\n",
    "## Preparing Small Molecule (ligand) for Docking\n",
    "Preparing the **small molecule** structure for the **docking procedure**. Converting the **PDB file** to a **PDBQT file** format (AutoDock PDBQT: Protein Data Bank, with Partial Charges (Q), & Atom Types (T)), needed by **AutoDock Vina**. <br><br>\n",
    "The process adds **partial charges** and **atom types** to every atom. Besides, the **ligand flexibility** is also defined in the information contained in the file. The concept of **\"torsion tree\"** is used to represent the **rigid and rotatable** pieces of the **ligand**. A rigid piece (**\"root\"**) is defined, with zero or more rotatable pieces (**\"branches\"**), hanging from the root, and defining the **rotatable bonds**.<br><br>\n",
    "More info about **PDBQT file format** can be found in the [AutoDock FAQ pages](http://autodock.scripps.edu/faqs-help/faq/what-is-the-format-of-a-pdbqt-file).\n",
    "\n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [babel_convert](https://biobb-chemistry.readthedocs.io/en/latest/babelm.html#module-babelm.babel_convert) from **biobb_chemistry.babelm.babel_convert**\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from biobb_chemistry.babelm.babel_convert import babel_convert\n",
    "\n",
    "prep_ligand = \"prep_ligand.pdbqt\"\n",
    "prop = {\n",
    "    \"input_format\": \"pdb\",\n",
    "    \"output_format\": \"pdbqt\",\n",
    "    \"obabel_path\": \"obabel\"\n",
    "}\n",
    "\n",
    "babel_convert(input_path = ligand,\n",
    "            output_path = prep_ligand,\n",
    "            properties=prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"viewDrug\"></a>\n",
    "### Visualizing small molecule (drug)\n",
    "Visualizing the desired **drug** to be docked to the **target protein**, using **NGL viewer**.<br>\n",
    "- **Left panel**: **PDB-formatted** file, with all hydrogen atoms.\n",
    "- **Right panel**: **PDBqt-formatted** file (AutoDock Vina-compatible), with **united atom model** (only polar hydrogens are placed in the structures to correctly type heavy atoms as hydrogen bond donors).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from ipywidgets import HBox\n",
    "\n",
    "v0 = nglview.show_structure_file(ligand)\n",
    "v1 = nglview.show_structure_file(prep_ligand)\n",
    "\n",
    "v0._set_size('500px', '')\n",
    "v1._set_size('500px', '')\n",
    "\n",
    "def on_change(change):\n",
    "    v1._set_camera_orientation(change['new'])\n",
    "    \n",
    "v0.observe(on_change, ['_camera_orientation'])\n",
    "\n",
    "HBox([v0, v1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='ngl5.png'></img>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"protein_pdb2pdbqt\"></a>\n",
    "***\n",
    "## Preparing Target Protein for Docking\n",
    "Preparing the **target protein** structure for the **docking procedure**. Converting the **PDB file** to a **PDBqt file**, needed by **AutoDock Vina**. Similarly to the previous step, the process adds **partial charges** and **atom types** to every target protein atom. In this case, however, we are not taking into account **receptor flexibility**, although **Autodock Vina** allows some limited flexibility of selected **receptor side chains** [(see the documentation)](https://autodock-vina.readthedocs.io/en/latest/docking_flexible.html).<br>\n",
    "\n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [str_check_add_hydrogens](https://biobb-structure-utils.readthedocs.io/en/latest/utils.html#utils-str-check-add-hydrogens-module) from **biobb_structure_utils.utils.str_check_add_hydrogens**\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_structure_utils.utils.str_check_add_hydrogens import str_check_add_hydrogens\n",
    "\n",
    "prep_receptor = \"prep_receptor.pdbqt\"\n",
    "prop = {\n",
    "    \"charges\": True,\n",
    "    \"mode\": \"auto\"\n",
    "}\n",
    "\n",
    "str_check_add_hydrogens(input_structure_path = pdb_protein,\n",
    "            output_structure_path = prep_receptor,\n",
    "            properties=prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"docking\"></a>\n",
    "***\n",
    "## Running the Docking\n",
    "Running the **docking process** with the prepared files:\n",
    "- **ligand**\n",
    "- **target protein**\n",
    "- **binding site box**<br>\n",
    "\n",
    "using **AutoDock Vina**. <br>\n",
    "\n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [autodock_vina_run](https://biobb-vs.readthedocs.io/en/latest/vina.html#module-vina.autodock_vina_run) from **biobb_vs.vina.autodock_vina_run**\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_vs.vina.autodock_vina_run import autodock_vina_run\n",
    "\n",
    "output_vina_pdbqt = \"output_vina.pdbqt\"\n",
    "output_vina_log = \"output_vina.log\"\n",
    "\n",
    "autodock_vina_run(input_ligand_pdbqt_path = prep_ligand,\n",
    "             input_receptor_pdbqt_path = prep_receptor,\n",
    "             input_box_path = output_box,\n",
    "             output_pdbqt_path = output_vina_pdbqt,\n",
    "             output_log_path = output_vina_log)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"viewDocking\"></a>\n",
    "### Visualizing docking output poses\n",
    "Visualizing the generated **docking poses** for the **ligand**, using **NGL viewer**. <br>\n",
    "- **Left panel**: **Docking poses** displayed with atoms coloured by **partial charges** and **licorice** representation.\n",
    "- **Right panel**: **Docking poses** displayed with atoms coloured by **element** and **ball-and-stick** representation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from ipywidgets import HBox\n",
    "\n",
    "models = 'all'\n",
    "\n",
    "v0 = nglview.show_structure_file(output_vina_pdbqt, default=False)\n",
    "v0.add_representation(repr_type='licorice', \n",
    "                        selection=models,\n",
    "                       colorScheme= 'partialCharge')\n",
    "v0.center()\n",
    "v1 = nglview.show_structure_file(output_vina_pdbqt, default=False)\n",
    "v1.add_representation(repr_type='ball+stick', \n",
    "                        selection=models)\n",
    "v1.center()\n",
    "\n",
    "v0._set_size('500px', '')\n",
    "v1._set_size('500px', '')\n",
    "\n",
    "def on_change(change):\n",
    "    v1._set_camera_orientation(change['new'])\n",
    "    \n",
    "v0.observe(on_change, ['_camera_orientation'])\n",
    "\n",
    "HBox([v0, v1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='ngl6.png'></img>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"selectPose\"></a>\n",
    "### Select Docking Pose\n",
    "Select a specific **docking pose** from the output list for **visual inspection**.\n",
    "<br>\n",
    "Choose a **docking pose** from the **DropDown list**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from Bio.PDB import PDBParser\n",
    "parser = PDBParser(QUIET = True)\n",
    "structure = parser.get_structure(\"protein\", output_vina_pdbqt)\n",
    "models = []\n",
    "for i, m in enumerate(structure):\n",
    "    models.append(('model' + str(i), i))\n",
    "    \n",
    "mdsel = ipywidgets.Dropdown(\n",
    "    options=models,\n",
    "    description='Sel. model:',\n",
    "    disabled=False,\n",
    ")\n",
    "display(mdsel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"extractPose\"></a>\n",
    "***\n",
    "## Extract a Docking Pose\n",
    "Extract a specific **docking pose** from the **docking** outputs. <br>\n",
    "\n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [extract_model_pdbqt](https://biobb-vs.readthedocs.io/en/latest/utils.html#module-utils.extract_model_pdbqt) from **biobb_vs.utils.extract_model_pdbqt**\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_vs.utils.extract_model_pdbqt import extract_model_pdbqt\n",
    "\n",
    "output_pdbqt_model = \"output_model.pdbqt\"\n",
    "prop = {\n",
    "    \"model\": mdsel.value + 1\n",
    "}\n",
    "\n",
    "extract_model_pdbqt(input_pdbqt_path = output_vina_pdbqt,\n",
    "             output_pdbqt_path = output_pdbqt_model,\n",
    "            properties=prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"pdbqt2pdb\"></a>\n",
    "***\n",
    "## Converting Ligand Pose to PDB format\n",
    "Converting **ligand pose** to **PDB format**. <br>\n",
    "\n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [babel_convert](https://biobb-chemistry.readthedocs.io/en/latest/babelm.html#module-babelm.babel_convert) from **biobb_chemistry.babelm.babel_convert**\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from biobb_chemistry.babelm.babel_convert import babel_convert\n",
    "\n",
    "output_pdb_model = \"output_model.pdb\"\n",
    "prop = {\n",
    "    \"input_format\": \"pdbqt\",\n",
    "    \"output_format\": \"pdb\",\n",
    "    \"obabel_path\": \"obabel\"\n",
    "}\n",
    "\n",
    "babel_convert(input_path = output_pdbqt_model,\n",
    "             output_path = output_pdb_model,\n",
    "            properties=prop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"catPdb\"></a>\n",
    "***\n",
    "## Superposing Ligand Pose to the Target Protein Structure\n",
    "Superposing **ligand pose** to the target **protein structure**, in order to see the **protein-ligand docking conformation**. <br><br>Building a new **PDB file** with both **target and ligand** (binding pose) structures. <br>\n",
    "\n",
    "***\n",
    "**Building Blocks** used:\n",
    " - [cat_pdb](https://biobb-structure-utils.readthedocs.io/en/latest/utils.html#module-utils.cat_pdb) from **biobb_structure_utils.utils.cat_pdb**\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from biobb_structure_utils.utils.cat_pdb import cat_pdb\n",
    "\n",
    "output_structure = \"output_structure.pdb\"\n",
    "\n",
    "cat_pdb(input_structure1 = pdb_protein,\n",
    "             input_structure2 = output_pdb_model,\n",
    "             output_structure_path = output_structure)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id=\"viewFinal\"></a>\n",
    "### Comparing final result with experimental structure \n",
    "Visualizing and comparing the generated **protein-ligand** complex with the original **protein-ligand conformation** (downloaded from the PDB database), using **NGL viewer**. <br>\n",
    "- **Licorice, element-colored** representation: **Experimental pose**.\n",
    "- **Licorice, green-colored** representation: **Docking pose**.\n",
    "<br>\n",
    "\n",
    "Note that outputs from **AutoDock Vina** don't contain all the atoms, as the program works with a **united-atom representation** (i.e. only polar hydrogens)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "view = nglview.NGLWidget()\n",
    "\n",
    "# v1 = Experimental Structure\n",
    "v1 = view.add_component(download_pdb)\n",
    "\n",
    "v1.clear()\n",
    "v1.add_representation(repr_type='licorice', \n",
    "                     selection='STI',\n",
    "                     radius=0.5)\n",
    "\n",
    "# v2 = Docking result\n",
    "v2 = view.add_component(output_structure)\n",
    "v2.clear()\n",
    "v2.add_representation(repr_type='cartoon', colorScheme = 'sstruc')\n",
    "v2.add_representation(repr_type='licorice', radius=0.5, color= 'green', selection='UNL')\n",
    "\n",
    "view._remote_call('setSize', target='Widget', args=['','600px'])\n",
    "view\n",
    "\n",
    "# align reference and output\n",
    "code = \"\"\"\n",
    "var stage = this.stage;\n",
    "var clist_len = stage.compList.length;\n",
    "var i = 0;\n",
    "var s = [];\n",
    "for(i = 0; i <= clist_len; i++){\n",
    "    if(stage.compList[i] != undefined && stage.compList[i].structure != undefined) {        \n",
    "       s.push(stage.compList[i])\n",
    "    }\n",
    "}\n",
    "NGL.superpose(s[0].structure, s[1].structure, true, \".CA\")\n",
    "s[ 0 ].updateRepresentations({ position: true })\n",
    "s[ 0 ].autoView()\n",
    "\"\"\"\n",
    "\n",
    "view._execute_js_code(code)\n",
    "\n",
    "view"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='ngl7.png'></img>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "<a id=\"questions\"></a>\n",
    "\n",
    "## Questions & Comments\n",
    "\n",
    "Questions, issues, suggestions and comments are really welcome!\n",
    "\n",
    "* GitHub issues:\n",
    "    * [https://github.com/bioexcel/biobb](https://github.com/bioexcel/biobb)\n",
    "\n",
    "* BioExcel forum:\n",
    "    * [https://ask.bioexcel.eu/c/BioExcel-Building-Blocks-library](https://ask.bioexcel.eu/c/BioExcel-Building-Blocks-library)\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
